%   multiloss_MC v.1 (03/12/2011)
%   ******************************************************************
%   This program creates the monte carlo results for Komunjer and Owyang
%       The program is designed to generate forecast error data using a recursive forecasting method
%       The forecast errors are then used to estimate the forecaster loss using both the separable and nonseparable
%       loss functions

%   ******************************************************************
%   this part of the program generates artificial data for komunjer and owyang
%   the actual data is generated from a trivariate ar(1)
%   the forecast's are generated from a loss function which minimizes the expected forecast error
%   the forecasts are generated recursively with an estimation period of 100


function [results,params,res] = multiloss_MC
warning off
%   User defined preferences
T0 = 550;   %   the total sample
t_f = 300;  %   the number of periods in the estimation sample
T_drop = 100;%   the number of data samples dropped off to remove initial conditions
%   ******************************************************************
global ERROR_SERIES INSTRUMENTS P;
options = optimset('Display','off'); 
%   parameters used to generate the actual data
P = 2;
p = P;
a = 0.5 * eye(3);
sig = eye(3);
alpha = [0.5;0.5;0.5];
c = [1;2;3];
for i = 1:2
    for j = 1:3
        for k = 1:2
            disp(['Create the parameters ',' i = ',num2str(i),' j = ',num2str(j),' k = ',num2str(k)])
            params(i,j,k).a = a;
            params(i,j,k).a(1,:) = params(i,j,k).a(1,:) + (i-1) * [0, 0.09, -0.06];
            params(i,j,k).a(2,:) = params(i,j,k).a(2,:) + (i-1) * [-0.12, 0, 0.4];
            params(i,j,k).a(3,:) = params(i,j,k).a(3,:) + (i-1) * [0.03, 0.4, 0];
            params(i,j,k).sig = sig;
            params(i,j,k).sig(1,2) = params(i,j,k).sig(1,2) + 0.03 * 3 * (k-1);
            params(i,j,k).sig(2,1) = params(i,j,k).sig(2,1) + 0.03 * 3 * (k-1);
            params(i,j,k).sig(1,3) = params(i,j,k).sig(1,3) - 0.01 * 3 * (k-1);
            params(i,j,k).sig(3,1) = params(i,j,k).sig(3,1) - 0.01 * 3 * (k-1);
            params(i,j,k).sig(2,3) = params(i,j,k).sig(2,3) + 0.9 * (k-1);
            params(i,j,k).sig(3,2) = params(i,j,k).sig(3,2) + 0.9 * (k-1);
            %   parameters for the loss function
            params(i,j,k).alpha = alpha;
            params(i,j,k).alpha = params(i,j,k).alpha + (j-2) * [0.0; 0.3; 0.3];
            tau2(i,j,k) = 2 * params(i,j,k).alpha(2,1) - 1;
            tau3(i,j,k) = 2 * params(i,j,k).alpha(3,1) - 1;
            disp(['Generate the data',' i = ',num2str(i),' j = ',num2str(j)])
            for n = 1:100
                disp([' n = ',num2str(n)])
                %   ******************************************************************
                x0 = normrnd(0,1,3,1); % initial conditions for the data generator
                x = zeros(3,T0);
                for t = 1:T0
                    e_t = mvnrnd(zeros(1,3),params(i,j,k).sig)';
                    x(:,t) = c + params(i,j,k).a * x0 + e_t;
                    x0 = x(:,t);
                end
                %   ******************************************************************
                %   generate the forecast data
                %   the forecaster minimizes the sum of the (time separable) loss over the estimation period
                es = zeros(3,T0 - T_drop - t_f);
                inst = zeros(3,T0 - T_drop - t_f);
                for t = 1 + T_drop:T0 - t_f - 1
                    %   data for the estimation period
                    xdata = x(:,t:t + t_f);
                    xm1 = x(:,t - 1:t + t_f - 1);
                    xm1 = [ones(1,size(xm1,2));xm1];
                    A0 = zeros(size(params(i,j,k).alpha,1),size(params(i,j,k).alpha,1)+1);
                    AA = fminunc(@(A) mvloss(A,xdata,xm1,params(i,j,k).alpha,p),A0,options);
                    %   compute the forecast errors
                    es(:,t) = x(:,t + t_f + 1) - AA * [1;x(:,t + t_f)];
                    inst(:,t) = x(:,t + t_f);
                end
                comp = tau2(i,j,k) * es(2,:) +  tau3(i,j,k) * es(3,:);
                %   ******************************************************************
                %   estimation
                disp(['Estimation',' i = ',num2str(i),' j = ',num2str(j),' k = ',num2str(k),' n = ',num2str(n)])
                ERROR_SERIES = es(:,1 + T_drop:end)';
                INSTRUMENTS = [ones(size(inst,2)-T_drop,1) (inst(:,1 + T_drop:end))'];
                results(i,j,k,n) = gmm_multi3ik(ERROR_SERIES,INSTRUMENTS);
                res(i,j,k,n).bias_a1 = results(i,j,k).alfa_1 - params(i,j,k).alpha';
                res(i,j,k,n).bias_a2 = results(i,j,k).alfa_2 - params(i,j,k).alpha';
            end,
        end,
    end,
end,










